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A  two-dimensional  (2D)  EIS  simulation  approach  is  developed  by  solving  a  SOFC  unit  cell  model 
with  imposed  sinusoidal  voltage  perturbations  at  different  frequencies.  The  transient  SOFC  unit  cell 
model  describes  the  intricate  interdependency  among  the  ionic/electronic  conduction,  multi-component 
species  transport,  electrochemical  reaction  processes  and  electrode  microstructure  as  well  as  the  cou¬ 
pling  processes  of  mass,  energy,  momentum  transport  within  flow  channels.  The  model  calculates  the 
local  transient  response  and  impedance  spectra  as  a  function  of  channel  position.  The  effects  of  the  reac¬ 
tion  depletion,  product  accumulation  as  well  as  the  temperature  variation  along  the  flow  channels  on  the 
EIS  spectra  are  numerically  simulated  with  a  counter-flow  mode.  The  results  show  that  the  convection- 
diffusion  process  along  the  flow  channel  has  significant  effects  on  the  low  frequency  half  circle  of  the 
impedance  spectra.  The  temperature  oscillations  accumulate  along  the  flow  channels,  and  then  affect  the 
current  responses  which  probably  lead  to  an  electro-thermal  impedance  effects. 

©  2012  Elsevier  B.V.  All  rights  reserved. 


1.  Introduction 

Electrochemical  Impedance  Spectroscopy  (EIS)  is  a  widely  used 
characterization  method  to  distinguish  processes  occurring  on  dif¬ 
ferent  time-scales  and  gain  insights  into  limiting  factors  limiting  of 
the  solid  oxide  fuel  cell  performance. 

Several  researchers  have  developed  or  are  currently  devel¬ 
oping  SOFC  impedance  models.  Fleig  and  Maier  [1]  studied  the 
influence  of  imperfect  contact  between  porous  electrodes  on  the 
impedance  by  using  finite  elements  modeling.  Bieberle  and  Gauck- 
ler  [2]  developed  a  physicochemical  model  to  describe  the  detailed 
elementary  electrochemical  reactions  within  the  electrode.  To  sim¬ 
ulate  the  electrochemical  impedance  spectra,  the  models  were 
solved  directly  through  the  State  Space  Modeling  (SSM)  approach, 
which  is  widely  used  in  control  theory  for  solving  complex  differen¬ 
tial  equations.  Bessler  [3]  presented  a  new  computational  method 
for  calculating  the  impedance  based  on  detailed  physicochemical 
models.  The  method  was  based  on  numerical  simulations  of  the 
transient  behavior  of  the  reaction  system.  The  impedance  was  then 
calculated  in  the  time  domain  from  the  simulated  periodic  response 
of  the  system,  maintaining  its  full  nonlinear  response.  Zhu  and 
Kee  [4]  developed  a  time-accurate  model  to  analyze  EIS  spectra 
in  anode-supported  button  cells  with  internal  methane  reform¬ 
ing.  A  common  feature  of  all  these  modeling  approaches  is  the 
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assumption  of  in-plane  uniformity  which  could  simplify  the  system 
to  one  dimension  in  the  electrode  thickness  direction. 

However,  Schneider  et  al.  [5,6]  found  that  the  second  loop  of 
impedance  can  be  related  with  the  propagation  of  the  oscillation 
of  product  concentration  along  the  flow  channel,  and  they  experi¬ 
mentally  verified  the  existence  of  the  two-dimensional  effects,  and 
thus  the  current  response  in  upstream  parts  of  the  cell  affects  the 
impedance  response  in  downstream  regions. 

This  paper  contributes  a  two-dimensional  (2D)  EIS  simu¬ 
lation  approach  by  solving  a  SOFC  unit  cell  transient  model 
with  imposed  sinusoidal  voltage  perturbations  at  different  fre¬ 
quencies.  The  transient  SOFC  unit  cell  model  described  the 
intricate  interdependency  among  the  ionic  conduction,  electronic 
conduction,  multi-component  species  transport,  electrochem¬ 
ical  reaction  processes  and  electrode  microstructure  within 
Membrane-Electrode-Assembly  (MEA)  as  well  as  the  coupling  pro¬ 
cesses  of  mass,  energy,  momentum  transport  and  charge  transfer 
within  flow  channel  and  inter-connectors.  Then,  the  effects  of  the 
reaction  depletion,  product  accumulation  as  well  as  the  temper¬ 
ature  variation  along  the  flow  channels  on  the  EIS  spectra  was 
numerically  simulated  at  counter-flow  operation  mode  for  the 
SOFC  unit  cell. 

2.  Model  development 

2.1.  Model  assumptions  and  calculation  domain 

Numerical  models  of  planar  SOFC  were  developed  by  cou¬ 
pling  governing  equations  of  charge,  mass,  momentum,  energy 
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Nomenclature 

c  concentration  (mol  nrr3 ) 

C  specific  double  layer  capacitance  (F m-2 ) 

D  diffusion  coefficient  (m2  s-1 ) 

E  activation  energy  (J  mol-1 ) 

F  frequency  (Hz) 

/  current  (A) 

i  current  density  (A  nr  2 ) 

z'o  exchange  current  density  (A  nrr 2 ) 

Mf  molecular  weight  of  species  i  (kg  mol-1 ) 

P  pressure  (Pa) 

Q.  source  term  of  charge  balance  equations  (A  nrr3 ) 

R  gas  constant  (8.314 J  mol-1 1<-1 ) 

Ri  source  term  of  mass  balance  equations  (kg  nrr  3  s-1 ) 

STpb  active  surface  area  per  unit  volume  (m2  nrr3 ) 

Ltpb  active  length  per  unit  volume  (m  nrr3 ) 
t  time  (s) 

T  temperature  (I<) 

V  electric  potential  or  cell  voltage  (V) 

VitVj  diffusion  volume 

Wj  mass  fraction  of  species  i 

xk  molar  fraction  of  gas-phase  species  k 

Y  admittance  (S) 

Z  coordination  number 

Greek  letters 

a  transfer  coefficient 

P  tuning  parameter  ( nr  2 ) 

s  porosity 

r)  overpotential  (V) 

p  density  (kg nr3) 

o  conductivity  (S  nr1 ) 

r  tortuosity 

oo  angle  frequency  (rad  s-1 ) 

Superscripts 
eff  effective 

Subscripts 

ac  anode  chamber 

act  activation/active 

an  anode 

avg  average 

bulk  bulk  phase 

ca  cathode 

cc  cathode  chamber 

elec  electronic 

F  Faraday 

inter  interface 

ion  ionic 

Kn  Knudsen 

sp  support 


transport  equations.  The  continuum  micro  model  for  electrode  sta¬ 
tistical  property  description  has  been  involved.  The  geometry  of 
this  SOFC  unit  is  shown  in  Fig.  1. 

The  geometry  of  the  calculation  domain  is  based  on  the  typical 
planar  SOFC  design  with  length  1 00  mm,  air  and  fuel  channel  height 
(rib  thickness)  1  mm,  interconnect  thickness  2.5  mm  (1 .5  mm  with¬ 
out  rib  part).  The  simplification  of  the  two-dimensional  calculation 
domain  would  be  reasonable  in  some  of  the  planar  SOFC  design 
where  the  rib-channel  width  ratio  or  channel  height-width  is  rel¬ 
atively  small. 


The  main  assumptions  made  in  this  model  are  shown  as  the 
following: 

(1)  The  electrochemical  reactions  spatially  occurred  along  elec¬ 
trode  thickness  within  electrode.  The  reaction  active  sites  are 
assumed  to  be  uniformly  distributed  in  each  electrode  layer. 
The  two  conducting  phases  (electronic  and  ionic)  are  consid¬ 
ered  to  be  continuous  and  homogeneous  in  each  layer. 

(2)  The  ionic  and  electronic  charge  transport  processes  take 
place  in  PEN  (Positive  electrode-Electrolyte-Negative  electrode 
assembly)  and  interconnect.  The  electronic  transport  in  inter¬ 
connect  as  well  as  the  contact  and  additional  resistance  due  to 
interconnect  were  not  considered. 

(3)  H2-H20-N2  mixtures  and  air  were  chosen  as  the  fuel  and  oxi¬ 
dant  of  SOFC  individually.  All  gas  mixtures  are  considered  as 
ideal  gases. 

(4)  Convection  flux  is  neglected  in  the  porous  electrode  compared 
to  diffusion.  Pressure  gradients  in  the  porous  electrode  are  also 
neglected. 

(5)  The  outer  boundaries  of  the  unit  cell,  for  example,  the  surface 
of  inter-connector,  the  inlet  and  outlet  of  fuel  and  air,  etc.,  were 
considered  as  thermal  insulated. 

2.2.  Governing  equations 

The  continuum  micro-scale  PEN  sub-model  coupled  the 
processes  of  ion/electronic  charge  transport,  porous  electrode  dif¬ 
fusion  as  well  as  the  effects  of  electrode  micro-structures.  And  in 
the  channel  there  are  only  mass  conservation  equation  and  energy 
conservation  equation. 

2.2 A.  Charge  conservation 

It  should  be  noted  that  ions  only  exist  in  the  PEN,  and  elec¬ 
trons  exist  in  the  anode,  cathode  and  inter-connectors.  Governing 
equations  for  transient  charge  balance  were  shown  as  follows  [  7,8  ] : 

•  Ionic  charge  balance  at  the  cathode: 


Qn,ca^act,cad(  Vjon?ca  Vel,ca) 

Wt 


+  V'(-CTfon,caWion,ca)  =  Qion,c 


(C™  f  aneF(Vel  ca  -  Vion  c  -  Vref,ca) 

=  -'0,caSact,ca  eXP  ' 


"02 


RT 


-  exp  - 


(1  a)neF(Ve j  ca  17jon  ca  17ref  ca) 


RT 


(1) 


Electronic  charge  balance  at  the  cathode: 


Ql,caSact,ca9(Vel,ca  —  Vion  ,ca) 

Ft 

+  V  •  (— <fcaWel,ca)  =  Qel.ca  =  -Qion.c 


(2) 


Ionic  charge  balance  at  the  anode: 


Qn,anSact,an^(  Vion,an  ^el,an) 

di 

—  Qion,an  —  *0,anSact,an 


+  V'(-<n,anWion,an) 
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Cathode  layer  (LSM/ScSZ) 
Electrolyte  layer  (Dense  ScSZ  ) 
Anode  active  interlayer  (Ni/ScSZ) 
Anode  support  layer  (Ni/YSZ) 


Calculation  domain 


an. 


inter/cc 


^ca/cc 

dQ 


elec/ca 


dQ, 

dQ 

dQ 


elec/an  act 


'elec/an  act 


an_sp/ac 


dQ, 

dQ:, 


inter/ac 


inter/cc 


TPB 


(Inter  connector/air  channel) 
(Cathode/air  channel) 
(Electrolyte/cathode) 
(Electrolyte/anode  active  layer) 
(Electrolyte/anode  active  layer) 
(Anode  support  layer/fuel  channel) 
(Inter  connector/fuel  channel) 

(Inter  connector/air  channel) 


Fig.  1.  SOFC  unit  cell  geometry  and  calculation  domain. 


rTPB 

%_ 

rbulk 

% 


exp 


aneF(Ve  i 

,an  ^ion,an  Vref,  an) 

RT 


rTPB 

h2o 

rbulk 

H20 


x  exp  - 


(1  ^)^e^(^el,an  ^ion,an  ^ref,an) 

RT 


(3) 


electrode  and  the  implicit  relationship  between  the  species  mole 
fraction  and  molar  diffusion  fluxes  was  provided, 


,i  =  1,2,3,  ...,n 


(5) 


where  Deff  is  the  effective  molecular  diffusion  coefficients  and  Dfff . 

ij  kn,t 

is  the  effective  Knudsen  diffusion  coefficients.  This  formula  can  be 
•  Electronic  charge  balance  at  the  anode :  written  as  the  explicit  relationship  between  species  model  fractions 

and  molar  fluxes  using  matrix  notation, 


Qn,anSact,an9(^el,an  ^ion,an) 

di 

+  V  •  (— CTei  an Wei?an)  =  Qel,an  —  —  Qion,an  (4) 

where  t  is  time,  Cdl  is  the  specific  interface  double-layer  capacitance 
between  electronic  and  ionic  conductor  phases;  Q.  is  the  transfer 
current  source;  Ve\  and  Vion  denote  the  electric  potential  of  the  two 
conductor  phases  [9]. 


2.2.2.  Mass  conservation 

The  dusty-gas  model  is  used  for  describing  diffusion  processes 
within  porous  electrode.  The  dusty-gas  model  could  be  formulated 
as  follows  while  neglecting  the  pressure  gradient  in  the  porous 


(N,)  =  -c[B]_1(Vx,)  =  -c[D](Vx,)  = 


(6) 


where  the  elements  of  square  matrix  [B]  can  be  formulated  as, 


1 


Bv  = 


Kn,eff 


D 


■E 


*k 


,MS,eff  ’ 


J=l 


Xi 


£)MS,eff  ’ 


k= 1  ik 

j^i 


i,j  =  1,2,3,  ...,n 


(7) 


The  [D]  is  the  inverse  of  the  matrix  [B],  the  elements  of  [D]  can  be 
treated  as  the  effective  multi  component  diffusivity  according  to 
the  dimension  of  its  unit. 
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Table  1 

Governing  equations  for  mass  conservation. 


Domain  Governing  equations 


Anode  3(£an^nWH2  j  +  v  .  |  -panwH2  ^^(D^jVXj)  +  panwH2u 


j= i 


-itrans,  an  ^act,  an  ^H2 


3(£anPanWH  0)  „  /  \ 

- dt  2  +  V  ■  -AmWH2o  y  (Pn2QjVxf)  +  AmWH2oU 

-Itrans,anSact,an^H20 


2F 


Cathode  9(£ca^aW°2 }  +  v  j  -pcaw02  ^^(Dq2  j  Vx o2 )  +  pcaw02  u 


itrans,  ca  Sact,  ca  ^02 


Fuel 

channel 


3(PanWH2 

dt 


+  v  [  -panWH2  y^(PH2j, chanel  VXj)  +  PanWHz  U  )  =0 


9tH20)  +  V  ■  Pan Wh2o^^(5h2OJ, channel  VXj)  +  panWH2oU ^  =  0 


Air 

channel 


3(PcaWQ2 ) 

dt 


+  v  -PcaWo2  y^(Pp2  J, channel  VXo2  )  +  PcaWp2  U  ]  =0 


Then,  we  can  get, 

3(epwf) 


at 


+  V  -pWf^^DijVXj  +  pwj-u  =  K, 


(8) 


The  equations  for  electrode  mass  balances  can  be  summarized 
as  in  Table  1. 


Table  2 

Governing  equations  for  momentum  conservation. 


Domain 

Governing  equations 

Anode 

pft  +v  («Pw)  =  -«Vp  +  v  [£/imix(V5  +  (Vi))T)]-  0f*e2v 

Cathode 

p\ f  +  V  •  ( spvv )  -  -s Vp  +  V  •  [«/xmix(Vt>  +  (VP)T)]  -  Ossa. s2v 

Fuel  channel 

p|  +  V  ■  (piiv)  =  -Vp  +  V  •  [Mmix(VS+  (VD)t)] 

Air  channel 

pf  +  V  •  (pvv)  =  -Vp  +  V  •  [/WV5  +  (V5)T)] 

Table  3 

Governing  equations  for  energy  balance. 

Domain 

Governing  equations 

Solid  structures 

W)  +  V  ■  (-AeffVr)  =  Qtex 

Fuel  and  air  channels  +  V  ■  (-AVT  +  pCpTu)  =  0 

2.2.4.  Energy  balance 

The  heat  transport  processes  within  SOFC  unit  cell  are  usually 
complex  as  shown  in  following: 


•  MEA:  conductive  and  radiative  heat  transfer,  and  the  heat  sources 
include  the  ohmic  heat  due  to  ionic  and  electronic  conducting,  the 
reversible  heat  effects  of  entropy  change  in  electrochemical  reac¬ 
tions,  as  well  as  the  irreversible  heat  generation  due  to  activation 
polarizations. 

•  Interconnect:  conductive  and  radiative  heat  transfer. 

•  Gas  channels:  convective  heat  transfer. 


The  governing  equations  for  energy  balance  are  shown  in 
Table  3. 

As  and  Ag  are  thermal  conductivities  of  solid  phase  and  gas  phase 
materials  in  porous  electrodes.  Qheat  is  the  heat  source  term,  and 
can  be  calculated  as  following  equation  in  different  solid  structures, 


2.2.3.  Momentum  conservation 

The  general  formulated  non-isothermal  weakly  compressible 
Navier-Stokes  equation  with  continuity  equation  was  used  to 
describe  the  momentum  balance  in  the  gas  channels, 

T  2  1 

p(Vu  +  (Vu)  )-  -rjV  u  , 

^  +  V.(pu)  =  0  (9) 

where  u  is  the  velocity  vector,  p  is  pressure,  r\  is  the  dynamic  viscos¬ 
ity.  p  is  the  density.  When  temperature  variation  was  considered, 
the  density  will  be  a  function  of  the  temperature,  the  momentum 
equations,  continuity  equation  and  temperature  equation  form  a 
closely  coupled  system  of  equations. 

In  porous  sub-domains,  the  flow  variables  and  fluid  properties 
are  defined  at  any  point  inside  the  medium  by  means  of  averaging  of 
the  actual  variables  and  properties  over  a  certain  volume  surround¬ 
ing  the  point.  The  flow  velocities  are  defined  as  superficial  volume 
averages,  and  they  correspond  to  a  unit  volume  of  the  medium 
including  both  pores  and  matrix.  They  are  sometimes  called  Darcy 
velocities,  defined  as  volume  flow  rates  per  unit  cross  section  of  the 
medium.  The  presented  approach  eliminates  the  need  for  explicit 
consideration  of  interfaces.  Thus,  the  flow  can  be  modeled  by  using 
the  same  unknown  variables  for  the  entire  domain  including  free 
flow  sub-domains  and  porous  sub-domains.  The  distinction  is  made 
via  switching  on  and  off  certain  terms  in  the  governing  equations. 
The  conservation  equations  are  summarized  in  Table  2. 


p—  +  pu  -Vu  =  -Vp  +  V  • 
dt 


{0,  inter-connector 
Qohm,  electrolyte 

Qohm  +  Qrev  +  Qi,r.  porous  electrode 


(10) 


where  the  Qohm  is  the  ohmic  heat  and  can  be  calculated  as, 

(  i2 

— ,  electrolyte 

Qohm  =  <  f°n  j2  HD 

— —  +  — ,  porous  electrode 
V  ^ion  ^el 

where  i  is  local  electronic  or  ionic  current  density,  a  is  electric 
conductivity. 

Qrev  is  the  reversible  heat  effects  of  entropy  change  in  electro¬ 
chemical  reactions,  and  can  be  calculated  as, 


itrans,an5act,an^(5H20  +  2Se-  —  5h2  —  ^o2-  ) 
2 F 


, anode 


itrans,caSact,ca^(SQ2-  —  (l/2)5o2  —  25e- ),  Cathode 
2 F 


(12) 


where  S  is  the  molar  entropies  of  species.  Then  the  anode  entropy 
change  can  be  calculated  by  subtracting  cathode  entropy  change 
from  whole  the  entropy  change  of  electrochemical  reaction.  Qirr 
is  the  irreversible  heat  generation  due  to  activation  polarizations 
which  can  be  formulated  as, 


Qirr  —  I^Qell  —  I^Qionl 


(13) 
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For  the  gas  channels,  the  energy  balance  equations  can  be  for¬ 
mulated  by  considering  convective  and  conductive  heat  transfer 
processes, 

+  v  •  (— A.VT  +  pCpTu)  =  0  (14) 

where  Cp  is  the  specific  heat  capacity,  and  A  is  thermal  conductivi¬ 
ties  of  gas  in  gas  channels. 

The  radiative  heat  transfer  between  the  MEA  solid  structure 
and  interconnect  has  been  considered  by  introducing  a  heat  flux 
of  surface-to-surface  radiation  in  heat  flux  boundary  conditions  at 
electrode/channel  and  interconnect/channel  interfaces. 


2.4.  Model  solution  and  EIS  spectra  synthesis 

The  problem  was  solved  for  a  given  cell  voltage.  The  perturba¬ 
tion  voltage  is  chosen  as  10  mV  to  assure  the  EIS  linearity,  and  in 
fact  this  value  is  set  as  the  same  with  that  for  EIS  experimental  tests. 
The  outputs  of  the  model  are  the  distributions  of  current  density 
and  species  concentrations.  The  calculations  were  performed  using 
the  finite  element  commercial  software  COMSOL  MULTIPEISICS®. 
Since  large  differences  exist  in  length  scales  of  MEA  structure  and 
channels/interconnect,  the  asymmetry  mapping  structured  quadri¬ 
lateral  mesh  was  used  in  the  solution.  The  mesh  consists  of  4940 
elements  with  the  number  of  freedom  degree  20,036  which  could 


4 rad,  electrode-inter  — 


°B  ( ^electrode  ^inter  ) 


( ( 1  /^electrode  )  1 )  +  ( 1  /^electrode-inter  )  +  CAlectrode  /Anter  )( ( 1  /^inter  )  1 ) 


(15) 


4 rad,  inter-electrode  — 


"Bin '-n 


electrode  / 


((1/^inter)  1 )  +  ( 1 /^inter-electrode  )  +  (Anter /Alectrode)((l  /^electrode)  1 ) 


(16) 


where  £  is  surface  emissivity.  crB  is  the  Stefan-Boltzmann  constant 
(5.67e-8  Wm-2 1<-4)  F  is  the  radiative  view  factor  which  can  be 
calculated  according  to  the  SOFC  geometry. 


2.3.  Boundary  conditions 

In  order  to  solve  the  system  of  coupled  partial  differential  equa¬ 
tions  for  charge,  mass,  momentum  and  energy  balances,  the  name 
of  boundary  conditions  are  specified  in  Fig.  1.  First  of  all,  the  up  and 
bottom  boundaries  9*f2inter/cc  were  set  as  periodical  for  each  phys¬ 
ical  field.  The  detailed  settings  of  other  boundaries  were  listed  as 
following: 

Ionic  charge  balance.  At  the  electrode/electrolyte  interfaces 
3if2Ca/eiec*  ^^an.act/eiec*  the  continuity  boundary  condition  was 
applied  to  specify  that  the  normal  ionic  current  are  continuous 
across  the  interior  boundary.  The  electric-insulation  condition, 
specifies  no  current  flow  across  the  boundary,  was  applied  at  other 
boundaries. 

Electronic  charge  balance.  At  the  electrode/channel  interfaces 
3^an_sp/ac.  3^ca/cc  the  electric  potentials  Vcell,an  and  Vcell,ca  were 
specified,  respectively.  Elere,  Vcell  an  was  set  as  zero  and  Vcellca  is 
set  as  the  cell  working  voltage.  At  the  electrode/electrolyte  inter¬ 
faces  9*f2Ca/cc*  d£?an_act/ac  and  other  boundaries,  electric-insulation 
conditions  were  specified. 

Mass  balance.  The  mass  fraction  was  specified  at  the  inlet  of 
fuel  and  air  channels  while  the  convective  flux  boundary  condi¬ 
tion,  which  indicates  that  the  mass  transport  out  of  the  domain 
is  dominated  by  convection,  was  specified  at  outlets  of  gas  chan¬ 
nels.  At  the  electrode/electrolyte  interfaces  9-f2ca/elec,  3^an_act/eiec 
and  other  boundaries,  the  gas  insulation  conditions  were  specified. 

Momentum  balance.  The  velocity  was  specified  at  the  inlets  of 
fuel  and  air  channels  while  the  pressure  condition  was  specified 
at  the  outlets  of  the  channels.  Other  boundaries  were  set  as  walls 
with  no  slip  which  indicates  that  the  velocity  of  the  gas  at  the  wall 
is  zero. 

Energy  balance.  The  gas  temperature  was  given  at  the  inlet  of 
fuel  and  air  channels.  The  convective  heat  flux  boundary,  which 
means  n(-AVT)  =  0,  was  specified  at  the  outlets  of  gas  chan¬ 
nels.  Energy  transport  processes  of  gas  and  solid  are  coupled  at 
the  electrode/channel  interfaces  and  interconnect/channel  inter¬ 
faces  by  specifying  that  the  temperature  are  continuous  across  the 
boundaries,  and  the  radiation  heat  flux  boundary  was  also  speci¬ 
fied  at  these  interfaces.  The  other  boundaries  were  specified  as  heat 
insulation. 


make  sure  the  calculation  resolution  independent  with  the  mesh 
element  number.  We  have  placed  76  elements  along  axial  direction 
and  65  elements  along  electrode  thickness  direction.  For  electrode 
thickness  direction,  we  use  5  elements  in  electrolyte,  cathode  and 
anode  active  layer,  and  use  12  elements  in  anode  channel,  cathode 
channel.  For  anode  support  layer  and  current  collector  layer,  13 
elements  were  used.  The  default  stationary  nonlinear  solver  (uses 
damped  Newton  method  with  direct  linear  system  solver  UMF- 
PACK)  of  COMSOL  MULTIPEISICS®  was  used  for  a  certain  operating 
case,  while  the  parametric  nonlinear  solver  was  used  to  get  the  EIS 
spectra.  The  relative  tolerance  of  the  model  was  specified  as  le-5. 

The  average  current  density  was  calculated  as  following  equa¬ 
tion  in  electrolyte  layer  with  constant  coordinate  y: 


/ 


hon,  local ^ 


(17) 


where  iion, local  is  the  local  ionic  current  density  in  electrolyte  layer 
and  L  is  the  length  of  the  calculation  domain  along  x  direction. 

In  order  to  generate  a  complete  EIS  spectrum,  the  model  needs 
to  be  run  for  several  periods  just  like  the  experimental  measure¬ 
ments  of  the  EIS.  The  time  step  is  equal  to  one  twentieths  periods, 
and  the  time  span  is  equal  to  15  periods.  The  detailed  calcula¬ 
tion  method  can  be  found  in  our  previous  studies  [10].  While,  the 
two-dimensional  EIS  simulation  is  a  time-consuming  problem.  We 
calculate  the  current  response  at  different  frequencies  separately 
using  four  small  computation  workstation  (Intel  Xeon  3.2  GEIz,  8GB 
Memory).  The  total  calculation  time  for  a  whole  EIS  spectrum  is 
around  2  h.  In  order  to  calculate  the  impedance  of  the  fuel  cell, 
the  model  is  simulated  at  the  selected  voltage  for  10  s  to  achieve 
steady  state.  Then,  the  model  initial  condition  was  set  as  the  steady 
state,  and  the  model  is  solved  repeatedly  with  sinusoidal  signals 
superimposed  on  the  voltage. 


3.  Results  and  discussion 

3.1.  Model  parameters 

From  the  model  development  section  of  charge  balance  and 
mass  balance  equations,  the  effective  reaction  area,  TPB  area  and 
TPB  length  are  needed  in  the  calculation.  The  parameter  Stpb  can 
be  formulated  by  using  the  particle  coordination  number  in  binary 
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random  packing  of  spheres  together  with  percolation  theory  as 

[11,12] 

n  sin2  6»re2|/Vtotne|nioZelZioPelPio 

jtpb  = - £ -  (loj 

where  Z  is  the  mean  coordination  number,  reJ  is  the  mean  radius 
of  the  electronic  conductor  particle,  6  is  the  contact  angle  between 
the  electronic  and  ionic  conductors  particles,  Nt0 1  is  the  total  num¬ 
ber  of  particles  per  unit  volume,  nel  and  nio  are  the  fraction  number 
of  electronic  and  ionic  conductor  particles,  Zel  and  Zio  are  the  coor¬ 
dination  numbers  of  electron  and  ion  conductor  particles,  and  Pel 
and  Pio  are  the  whole  range  connection  probabilities  of  the  same 
kind  particles.  Here  Sjpb  is  treated  as  the  contact  area  of  different 
types  of  particles. 

Similar  to  the  calculation  Sjpb,  the  parameter  LTPB  (TPB  length 
per  unit  volume)  can  be  calculated  as, 

Ltpb  =  2n rel  sin  PNtotnelZel_ioPelPio  (19) 

The  pore  structure  of  the  anode  support  layer  was  charac¬ 
terized  using  mercury  porosimeter.  The  mean  pore  diameter, 
porosity  and  total  pore  area  were  found  to  be  0.387  |jim,  0.335 
and  8.54  x  106m2m-3,  respectively.  To  simplify  the  calculation, 
the  mean  particle  diameters  of  the  two  conductors  are  assumed 
to  be  the  same  and  equal  to  the  mean  pore  diameter  [13].  With  this 
assumption  and  the  expressions  of  effective  reaction  areas,  the  cal¬ 
culated  value  of  Sjpb  in  the  anode  support  layer  is  2.22  x  105.  The 
mean  pore  size  in  anode  support  layer  is  assumed  1 .2  and  1 .5  times 
of  that  in  cathode  layer  and  anode  active  layer  based  on  SEM  image 
from  quantitative  stereology. 

The  physical  properties  of  gas  mixture,  e.g.  specific  heat  capac¬ 
ity,  heat  conductivity,  viscosity,  et  al.  were  calculated  in  the 
well-correlated  equations  with  temperature  in  the  works  of  Young 
et  al.  [14].  The  physical  properties  of  solid  structures  are  listed  in 
Table  4.  These  parameter  values  were  based  on  the  values  from 
published  literatures  and  the  temperature  dependence  of  the  heat 
conductivity,  density  and  heat  capacity  of  all  solid  materials  were 
all  considered  negligible  in  the  temperature  range  concerned  here. 


Table  4 

Model  parameters. 


Parameters 

Value 

Electric  conductivity  (S  m- 1 ) 

Ionic  conductor  ScSZ,  crScsz 

6.92e4  exp(  -9681/7) 

Ionic  conductor  YSZ,  crYsz 

3.34e4exp(-10, 300/7) 

Electronic  conductor  LSM,  <tLsm 

4.2e7/7exp(-l  150/7) 

Electronic  conductor  Ni,  <rNi 

3.27e6-1065.37 

Equivalent  ionic  conductivity  of  electrolyte 

0.0027-1.4483 

layer 

Heat  conductivity  (W  m- 1  K-1 ) 

Anode  (support/active  layers)  A,an_Sp,  A,an_act 

6.23 

Cathode,  A.ca 

9.6 

Electrolyte,  A.eiec 

2.7 

Interconnect,  A.inter 

25 

Density  (kgm-3) 

Anode  (support/active  layers),  pan_Sp,  pan_ act 

6870 

Cathode,  pca 

6570 

Electrolyte,  peiec 

2000 

Interconnect,  p^ter 

8000 

Heat  capacity  (J  kg-1  K_1 ) 

Anode  (support/active  layers),  Cp,an_sp,  Cp,an.act 

420 

Cathode,  Cp>ca 

390 

Electrolyte,  CPieiec 

300 

Interconnect,  Center 

500 

Emissivity 

Anode 

0.8 

Cathode 

0.8 

Interconnect 

0.1 

Parameters  of  electrochemical  reaction  kinetics 

Interface  conductivity  Am_sP/^an_act/&a 

6.8el2/6.8el2/5.8el0 

(C2_1  nrr2) 

Activation  energy  Fact, an_sp /fact, an_act/Fact,ca 

120,000/120,000/130,000 

(J  mol-1 ) 

Microstructures 

Anode  support  layer 

0.364/16.5/0.48 

porosity/tortuosity/pore  radius  (|jim) 

Anode  active  layer  porosity/tortuosity/pore 

0.364/16.5/0.31 

radius  (p,m) 

Cathode  porosity/tortuosity/pore  radius 

0.364/3/0.4 

(pan) 

3.3.  Two  dimensional  EIS  spectra  modeling  at  constant 
temperature 


3.2.  PEN  sub  model  validation 

The  button  cell  experimental  data  was  used  to  validate  the  PEN 
sub-model  by  neglecting  the  effect  of  momentum  and  heat  transfer. 
The  test  setup,  validation  results  can  be  found  in  our  previous  paper 
[7-10].  The  parameters,  include  charge  transfer  reaction  kinetic 
parameter,  tortuosity  and  interface  capacity  for  both  cathode  and 
anode  were  well  validated  using  experimental  polarization  curves 
obtained  on  a  SOFC  button  cell  and  a  symmetric  cathode  cell,  whose 
material  and  structure  are  exactly  kept  the  same  with  the  MEA 
in  the  unit  cell  in  this  paper.  The  experimental  data  used  in  the 
model  validation  include  the  I-V curves  and  EIS  spectra  at  following 
operating  conditions: 

•  Different  cell  temperatures  at  750  °C,  800  °C  and  850  °C  with  oxy¬ 
gen  or  air  oxidant 

•  N2  dilution  at  800  °C  and  850  °C  with  N2  content  in  the  dry  fuel 
mixture:  0,  20%,  40%,  60%  and  80%. 

Fig.  2  shows  the  comparison  results  between  simulated  and 
experimental  polarization  curves  and  EIS  spectra.  It  can  be  seen 
that  the  PEN  EIS  model  can  well  predict  the  experimental  data  at 
different  temperatures  and  with  different  fuel  compositions.  These 
validations  will  be  the  basis  of  the  unit  cell  two-dimensional  EIS 
model. 


3.3  A.  Effects  of  operating  voltage  on  2D  EIS  spectra 

For  the  convenience  of  comparison  between  cases  with  different 
working  status,  one  set  of  operating  conditions  of  cell  referred  to 
“base  case”  was  selected  as  shown  in  Table  5, 

Fig.  3  shows  the  calculated  parameter  distributions  of  SOFC  unit 
cell  cross  section  at  0.7  V  and  at  constant  temperature  800  °C  in  the 
base  case.  It  can  be  seen  in  Fig.  3(a)  that  the  ionic  current  density 
decreases  along  flow  channel  direction.  The  current  density  value  at 
the  position  of  fuel  inlet  is  around  5000  A  m-2  which  is  near  the  cor¬ 
responding  button  cell  performance  as  shown  in  Fig.  2,  and  the  local 
current  density  gradually  decreases  to  around  1000 Am-2  at  the 
anode  outlet.  As  shown  in  Fig.  3(b)  the  H2  molar  fraction  decrease 
fast  along  flow  channel  due  to  the  continuous  H2  electrochemical 
consumption,  and  thus  the  local  cell  reversible  voltage  decease  and 
the  local  activation  and  concentration  polarization  losses  increase. 


Table  5 

Model  parameters  at  base  case. 


Operating  conditions 


Anode  Cathode 


Gas  inlet  mole  fraction 

Gas  inlet  mole  flowrate  (mol  h-1 ) 

Gas  inlet  velocity  (m  s-1 ) 

Gas  outlet  pressure  (Pa) 

Flow  configuration 
Cell  voltage 


96%  H2,  4%  H20  21%  02,  79%N2 

0.0225  0.58 

0.008  0.208 

101,325  101,325 

Counter  flow 
0.7 
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Z  1(0  cm2)  Current  density/(Am*2) 

Fig.  2.  Experimental  and  simulated  results  comparison:  (a)  EIS  spectra  at  750/800/850  °C;  (b)  polarization  curves  at  750/800/850  °C;  (c)  EIS  spectra  with  different  fuel 
compositions:  (d)  polarization  curves  with  different  fuel  compositions. 


Since  the  air  flow  rate  is  far  larger  than  the  fuel  flow  rate,  the  02 
molar  fraction  varies  slightly. 

Fig.  4  shows  the  calculated  local  impedance  spectra  in  the  base 
case  at  different  operating  cell  voltages  at  constant  temperature 
800  °C.  X  indicated  dimensionless  x-axis  along  fuel  flow  direc¬ 
tion.  The  impedance  increases  from  anode  gas  inlet  towards  anode 
gas  outlet  due  to  the  increased  charge  transfer  and  mass  trans¬ 
fer  impedance  due  to  fuel  consumption.  It  can  be  deduced  that 
the  average  impedance  can  predict  the  two  impedance  half  circles 
which  indicated  two  types  of  limiting  processes,  but  it  is  hard  to 
represent  the  impedance  status  at  different  positions  along  flow 
channel.  Especially  at  the  downstream  of  the  flow  channel,  the 
low  frequency  impedance  half  circle  increase  fast  due  to  the  mass 
transfer  in  both  electrode  thickness  direction  and  in  flow  channel 
direction  affected  by  the  product  accumulation.  These  results  indi¬ 
cated  that  ID  impedance  model  is  not  enough  for  using  in  unit 
cell  impedance  spectra  interpretation,  2D  effect,  e.g.  convective- 
diffusion  impedance,  has  to  be  considered.  Also,  it  can  be  seen  that 
the  radius  of  impedance  half  circle  shrinks  from  0.9  V  to  0.7  V  when 
X  is  larger  than  0.75  due  to  the  decrease  of  activation  polarization 
resistance.  At  this  time,  the  effect  of  diffusion  limitation  is  not  obvi¬ 
ous.  While  for  X  =  1,  the  low  frequency  impedance  for  cell  voltage 
at  0.7  V  increase  fast  due  to  the  increase  of  concentration  losses  at 
the  flow  channel  down  stream. 


3.3.2.  Effects  of  fuel  flow  rate  on  2D  EIS  spectra 

Fig.  5  shows  the  calculated  the  local  impedance  spectra  at  con¬ 
stant  temperature  800  °C  with  cell  voltage  at  0.7  V  by  increasing  the 
fuel  flow  rate  for  4  times  as  in  the  base  case  from  0.0225  molh-1 
to  0.09  mol  h-1.  The  results  indicated  that  the  low  frequency 
impedance  half  circle  is  very  small  at  fuel  inlet,  and  increase  gradu¬ 
ally  along  flow  direction.  In  this  case,  the  fuel  utilization  is  only  kept 
at  0.22  due  to  the  relatively  large  fuel  flow  rate  and  this  means  that 


the  fuel  is  not  consumed  much.  Thus,  it  can  be  seen  that  the  low 
frequency  impedance  half  circle  does  not  have  significant  increase 
as  shown  in  the  Fig.  4.  This  suggests  that  the  average  impedance 
can  be  good  prediction  in  the  unit  cell  impedance  only  when  the 
fuel  and  air  utilization  are  small. 


3.4.  Two  dimensional  EIS  spectra  modeling  by  coupling 
electro -thermal  effects 

3.4.1.  Electro -thermal  impedance  effects 

When  the  EIS  spectra  tests  are  carried  out  on  the  button  cell,  the 
operating  condition  is  kept  at  well  constant  temperature  oven.  Thus 
the  current  variation  in  EIS  tests  is  small,  and  the  effects  of  temper¬ 
ature  variation  due  to  the  cell  current  variation  can  be  neglected. 

However,  if  considering  a  SOFC  unit  cell  or  a  SOFC  stack,  the 
temperature  effects  from  AC  current  can  possibly  be  significant. 
Thus,  we  further  investigate  the  local  impedance  spectra  in  differ¬ 
ent  positions  along  flow  channel  with  consideration  of  temperature 
response  as  shown  in  Fig.  6.  The  results  indicated  that  the  heat 
transfer  will  have  significant  effects  on  the  impedance  spectra 
especially  at  low  frequency  part.  From  local  impedance  spectra  cal¬ 
culation,  the  inductive  half  circle  appeared  in  the  low  frequency 
part.  The  effect  of  temperature  variation  is  not  significant  near 
the  fuel  inlet.  But  the  effects  gradually  increase.  Before  X=  0.5,  the 
effects  are  mainly  represented  as  inductive  loops.  When  at  the 
down  stream  of  fuel  channel  (e.g.  X>0.5),  we  could  see  the  third 
impedance  half  circle  in  the  spectra.  The  impedance  leaded  by 
electro-thermal  effects  increase  to  the  value  compatible  to  the  low 
frequency  impedance  half  circle  (mass  transfer  impedance). 

To  further  understand  the  effects  of  temperature  variation  on 
EIS  spectra,  the  qualitative  analysis  of  impedance  characteristics 
can  be  done  when  considering  state  variables  include  the  operating 
voltage  V,  reacting  gas  concentration  C  and  temperature  T.  In  fact, 
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Current  density 
(A  m  2) 


0.02  0.04  0.06  0.08 

Position  along  flow  channel  direction(m) 


0  0.02  0.04  0.06  0.08  0.1 

Position  along  flow  channel  direction(m) 


Fig.  3.  Calculated  parameter  distributions:  (a)  ionic  current  density:  (b)  H2  molar 
fraction:  (c)  02  molar  fraction  at  0.7  V  and  at  constant  temperature  800 °C  in  the 
base  case. 


when  at  constant  operating  pressure,  by  assuming  at  the  steady 
sate,  the  Faraday  current  Jp  can  be  written  as  [15]: 

IF  =f(V,  T,  Q)  i=  1,2,3,...  (20) 

When  the  EIS  spectra  is  located  near  the  fuel  inlet  (e.g.  non- 
dimensional  position  Z=  0  along  flow  channel  direction).  Since 
the  fuel  concentration  is  high  at  fuel  inlet,  in  order  to  simplify 


Fig.  4.  Calculated  local  impedance  spectra  at  0.7  V  and  0.9  V  and  at  constant  tem¬ 
perature  800  °C  (1  mHz  to  1  MHz). 


Fig.  5.  Calculated  local  impedance  spectra  at  0.7  V  with  increasing  fuel  flow  rate  at 
constant  temperature  800  °C  (1  mHz  to  1  MHz). 


Fig.  6.  Calculated  global  and  local  impedance  spectra  in  the  base  case  at  0.7  V  con¬ 
sidered  energy  balance  (1  mHz  to  1  MHz). 
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the  discussions,  the  state  variables  reacting  gas  concentration  C 
could  be  assumed  to  be  constants.  Thus,  the  Faraday  current  is  con¬ 
sidered  to  the  function  of  state  variables  V  and  T.  Then,  in  order  to 
get  the  expression  of  Faraday  admittance  YF,  the  JF  can  be  Taylor 
expanded  to  linear  order  as: 


A IF 


AV  + 

T 


(21) 


Next,  to  express  A IF  as  the  only  function  of  AV,  the  AT  needs  to  be 
expressed  as  the  function  of  AV, 


A 


( 


dr 

at 


) 


=  JcdAT  = 


AT  ((d/dV)(dT/dt))v 
AV  jco-{{d/dT){dT/dt))v 


Position  along  flow  channel  direction(m) 


Then,  the  faradic  admittance  YF  can  be  obtained  as, 


y  =  A/f  / dl¥\  f  9/F\  AT  1  B 

F  av  ^av;T+  ^aT;vAv  *  j©+a 

where  the  parameters  A,  B  and  Rt  are  formulated  as: 


(24) 


(25) 

(26) 
(27) 


To  satisfy  the  stability  constraint  of  the  electrochemical  system, 
the  value  of  (3/3T)(3T/3t)v  should  be  negative,  which  means  the 
parameter  A  should  be  positive.  Rt  is  also  a  positive  parameter  since 
the  V  is  always  increase  with  JF  at  certain  temperature,  and  this 
represent  the  process  within  electric  double  layer.  And,  the  sign  of 
B  will  depend  on  the  values  of  m  and  b. 

It  is  easy  to  judge  that  the  value  of  m,  represented  as  (3/F/3T)v, 
should  be  positive  since  the  Faraday  current  increase  with  the  tem¬ 
perature  at  certain  operating  voltage.  While,  it  is  not  easy  for  the 
value  of  b  directly.  Since  the  electrochemical  system  is  operated  at 
steady  state,  thus  (3 T/3t)j  should  be  equal  to  zero.  Then,  we  have, 


dT  _  _dV  /3T\ 
dt~~dt[WJT 

Then,  this  formula  can  be  transformed  as: 

dT  _  ((3/3V)(3T/3t))r  _  b 

dV  ~  ~((3/3T)(3T/3t)V  “  A 


(28) 


(29) 


b-A%  (30) 

When  V  increases,  the  current  and  the  temperature  increase.  Thus, 
b  will  have  the  same  sign  with  A,  and  is  positive  in  this  situation. 
And,  this  also  indicates  that  the  B>  0.  Then,  the  admittance  YF  can 
be  further  formulated  as: 


Rt  +  ja)  +  A  Rt+  jcoL  +  R0  ^ 

where  R0  =A/B  >  0,  L  =  1  /B  >  0.  In  fact,  this  admittance  is  the  expres¬ 
sion  of  a  equivalent  circuit  with  the  parallel  connection  of  a  resistor 
Rt  and  a  complex  element  with  serious  connection  of  R0  and  L.  The 
circuit  formula  should  be  ( Rt{R0L )).  Thus,  according  to  above  dis¬ 
cussion,  it  is  obvious  that  when  considering  thermal  effects,  the 
inductive  loop  will  be  observed  in  the  EIS  spectra. 


Fig.  7.  Calculated  temperature  distribution  considering  thermal  effects  at  steady 
state  in  the  base  case. 


When  considering  the  local  EIS  spectra  near  the  fuel  outlet,  the 
concentration  polarization  will  became  significant  and  the  analyti¬ 
cal  analysis  will  be  more  difficult  [15].  While,  in  order  to  simplify  the 
analysis,  if  we  consider  an  extreme  condition  in  which  the  effects  of 
gas  concentration  is  far  significant  than  that  of  temperature  at  the 
downstream  of  flow  channel,  the  impedance  will  totally  represent 
the  large  low  frequency  capacitive  reactance  circle.  Thus,  we  could 
further  deduced  that  the  whether  the  local  EIS  spectra  is  capac¬ 
itance  circle  or  inductance  will  depend  on  the  complex  coupling 
interaction  among  charge  transfer  reaction,  gas  transport  process, 
temperature  variation.  And,  in  fact,  this  is  one  of  the  important  rea¬ 
son  to  develop  a  2D  mechanistic  EIS  model  in  present  manuscript. 
And,  from  above  preliminary  calculation  and  theoretical  discus¬ 
sions,  it  should  be  noted  that  the  SOFC  unit  cell  impedance  contains 
the  effects  which  is  hard  to  be  interpreted  by  common  zero-  or 
one-dimensional  EIS  spectra  models. 


3.4.2.  Temperature  distribution  and  response 

Fig.  7  shows  the  calculated  temperature  distribution  consider¬ 
ing  thermal  effects  at  steady  state  in  the  base  case.  The  fuel  and 
oxidant  temperature  increase  along  flow  channel  near  fuel  inlet, 
and  the  maximum  temperature  is  located  at  around  1/3  flow  chan¬ 
nel  length  from  fuel  inlet.  The  effects  from  temperature  variation 
may  lead  to  misinterpretation  of  the  diffusive  losses.  In  this  sit¬ 
uation,  2D  EIS  spectra  model  will  be  significantly  important.  In 
addition,  it  suggests  that  researchers  should  pay  more  attention 
in  unit  cell  impedance  experimental  measurements  to  avoid  the 
electro-thermal  coupling  effects,  and  to  clarify  the  contributions  to 
the  impedance  from  different  reaction  and  transport  processes. 

Fig.  8  shows  the  calculated  current  density  and  temperature 
responses  in  base  case  operating  condition.  It  can  be  seen  that  in 
this  case  the  temperature  is  not  a  constant  but  shows  a  periodic-like 
response  like  current  density.  In  fact,  the  temperature  and  current 
density  response  shows  the  coupling  effects  from  temperature  vari¬ 
ation  and  accumulation  effects  at  flow  channel  downstream.  When 
the  Z=  0,  it  can  be  seen  that  the  temperature  response  is  relatively 
smooth  since  only  the  energy  balance  is  coupled.  While,  when  the 
Z=0.5  and  Z=  1,  the  temperature  shows  irregular  response  due 
to  the  accumulation  effects  along  the  flow  channel.  Elowever,  the 
response  of  both  temperature  and  current  density  are  still  periodic. 
This  phenomenon  suggests  that  the  temperature  variation  affect 
the  impedance  spectra  shapes. 
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Dimensionless  position  along  flow  channel  direction:  Z  =  0 
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Dimensionless  position  along  flow  channel  direction:  Z  =  0.5 
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Fig.  8.  Calculated  current  density  and  temperature  response  in  the  base  case  at 
different  position  along  flow  channel  direction:  (a)  Z  =  0;  (b)  Z= 0.5  and  (c)  Z=  1 . 


3.4.3.  Effects  of  operating  voltage  on  electro-thermal  impedance 
Fig.  9  further  predicts  the  effects  of  operating  voltage  on  electro¬ 
thermal  impedance.  The  results  indicate  that  the  third  impedance 
half  circle  due  to  thermal  effects  at  0.7  V  is  larger  than  that  at  0.9  V. 
In  fact,  at  lower  cell  voltage,  the  current  density  is  relatively  larger, 
and  thus  the  thermal  effects  from  polarization  heat  became  signif¬ 
icant.  Especially  at  the  downstream  of  the  flow  channel,  the  third 
impedance  half  circle  increase  fast  due  to  the  coupling  temperature 


Fig.  9.  Effects  of  operating  cell  voltage  on  local  electro-thermal  impedance  spectra 
(1  mHz  to  1  MHz). 


variation  and  mass  transfer  in  both  electrode  thickness  direction 
and  in  flow  channel  direction  affected  by  the  product  accumula¬ 
tion.  These  results  further  indicated  that  the  periodical  temperature 
swing  due  to  the  perturbation  voltage  and  current  density  will  lead 
to  unusual  electro-thermal  impedance  effects.  While,  this  effect 
will  became  significant  at  relatively  larger  current  density  opera¬ 
tion  mode,  and  it  is  very  hard  to  be  interpreted  by  an  circuit  element 
in  a  equivalent  circuit  model  or  ID  impedance  model.  These  fur¬ 
ther  demonstrate  the  priority  of  direct  electrochemical  impedance 
modeling  approach  by  integrating  the  effects  of  multi-dimensional 
cell  geometry  and  multi-physics  coupling. 


4.  Conclusions 

In  this  paper,  a  two-dimensional  (2D)  EIS  simulation  approach 
was  developed  by  solving  a  SOFC  unit  cell  transient  model 
with  imposed  sinusoidal  voltage  perturbations  at  different  fre¬ 
quencies.  The  transient  SOFC  unit  cell  model  described  the 
intricate  interdependency  among  the  ionic  conduction,  electronic 
conduction,  multi-component  species  transport,  electrochem¬ 
ical  reaction  processes  and  electrode  microstructure  within 
Membrane-Electrode-Assembly  (MEA)  as  well  as  the  coupling  pro¬ 
cesses  of  mass,  energy,  momentum  transport  and  charge  transfer 
within  flow  channel  and  inter-connectors.  The  main  conclusions 
from  2D  EIS  simulation  results  include: 

The  convection-diffusion  process  along  the  flow  channel  has  sig¬ 
nificant  effects  on  the  low  frequency  half  circle  of  the  impedance 
spectra.  The  global  impedance  is  hard  to  represent  the  impedance 
status  at  different  positions  along  flow  channel. 

The  temperature  oscillations  accumulate  along  the  flow  chan¬ 
nels,  and  then  affect  the  current  response  which  probably  leads  to 
electro-thermal  impedance  effects. 

Equivalent  circuit  or  1 D  EIS  models  are  not  sufficient  to  describe 
the  SOFC  unit  cell  impedance  response,  and  the  work  in  this  study 
can  be  helpful  in  interpreting  EIS  spectra,  as  illustrated  in  identify¬ 
ing  the  contribution  of  different  reaction  and  transport  limitations. 
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